Heat wave propagation in a nonlinear chain 



0^ 
O 
O 
(N 

o 



O 
> 

ON 

o 

m 



00 

o 



X 



Francesco Piazza^ 'B and Stefano Lepri^'[l| 

^ Laboratoire de Biophysique Statistique, Institut de theorie des phenomenes physiques, 
Ecole Polytechnique Federale de Lausanne (EPFL) BSP-722, CH-1015 Lausanne, Switzerland 
^Istituto dei Sistemi Complessi, Consiglio Nazionale delle Ricerche, 
via Madonna del Piano 10, L50019 Sesto Fiorentino, Italy 
(Dated: January 9, 2009) 

We investigate the propagation of temperature perturbations in an array of coupled nonlinear 
oscillators at finite temperature. We evaluate the response function at equilibrium and show how 
the memory effects affect the diffusion properties. A comparison with nonequilibrium simulations 
reveals that the telegraph equation provides a reliable interpretative paradigm for describing quan- 
titatively the propagation of a heat pulse at the macroscopic level. The results could be of help in 
understanding and modeling energy transport in individual nanotubes. 
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I. INTRODUCTION 

The propagation of temperature fluctuations under 
nonequilibrium conditions may display significant devi- 
ations from a simple diffusive behavior. As it is know ^, 
heat can propagate in a wavelike form on time scales 
comparable with some typical relaxation time. This 
is accounted for by suitable extensions of the standard 
Fourier's law to include, for example, a finite response 
time for the heat current. On the other hand, it is of 
primary interest to derive such a phenomenological laws 
from microscopic dynamics. In this respect, much work 
has been performed recently, mostly on simplified mod- 
els, with the goal of understanding nonequilibrium heat- 
conducting states from first principles — . Much emphasis 
has been put on violations of Fourier's law in the station- 
ary case — as well as in describing anomalous heat diffu- 
sion — . Transient heat propagation has received much 
less attention—!^. As a matter of fact, the spreading of 
the energy perturbation field, yield complementary in- 
formation on how heat propagates through the systemS, 
and provides useful insight on the nature of heat carrying 
excitations. 

Besides those fundamental issues, there is a growing 
interest in understanding how heat is transported at the 
nanoscale — . In this context, carbon nanotube materi- 
als are of special relevance due to their exceptionally 
high thermal conductivity related to their quasi one- 
dimensional vibrational structure — i^. Recently, Osman 
and Srivastava— investigated the propagation of intense 
heat pulses in single-walled carbon nanotubes by means 
of molecular dynamics. They observed that, together 
with pulses traveling with the sound velocity, a secondary 
and slower peak can also propagate as a " second sound" 
type of wave — . Shiomi and Maruyamai^.^ argued that 
this should be related to relatively fast optical phonons, 
that due to the quasi one-dimensional structure, have 
very long relaxation times and thus contribute to wave- 
like heat propagation. 

It is thus relevant to study simplified models that can 
help to understand better the conduction properties. In 



this respect, an example is the length dependence of con- 
ductivity in carbon and boron-nitride nanotubes which 
has been very recently observed experimentally^'^. In- 
deed, experimental data are compatible with scaling laws 
theoretically predicted for simple nonlinear models. In 
this spirit, in the present paper we analyze the heat pulse 
propagation in a classical lattice model, a chain of classi- 
cal, harmonically coupled, oscillators with a quartic pin- 
ning (on-site) potential. In particular, we compare the 
linear-response predictions with nonequilibrium molecu- 
lar dynamics. 

The paper is organized as follows. In Sec. |TT]we intro- 
duce the microscopic one-dimensional model. After some 
general considerations about linear response we present 
the response functions as computed from molecular dy- 
namics simulations (Sec. IIII[) . Based on the numerical 
data, we introduce an approximated (single-pole) form 
for the response that leads to the telegraph equation for 
the temperature-field evolution. This equation, along 
the phenomenological values of its parameters, is pre- 
cious to compare with nonequilibrium simulations. The 
transient evolution of a temperature pulse is investigated 
in Sec. |TVl Qualitative and quantitative deviations from 
the behavior expected for the simple telegraph equation 
are highlighted. 



II. THE NONLINEAR CHAIN 

The model consists of an anharmonic chain of N par- 
ticles (each with unit mass) whose displacements are de- 
noted by Ui'. 



(1) 



In the linear approximation, where the cubic force term is 
dropped, the eigenfrequencies 17 of the associated normal 
modes are expressed as a function of the wavenumber q 

by 



fl^iq) = l + 2C(l-cosg) 



(2) 
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The units have been fixed in such a way to have a uni- 
tary gap in the spectrum. With this choice the only 
free parameters of the model are the energy per particle 
e = E/N and the coupling constant C. 

Several numerical studies (see e.g. Ref. [13 and [H as 
well as Ref. d and the bibliography within) clarified, 
that for models like (H]) the thermal conductivity is finite, 
and that a diffusive heat propagation is expected at long 
times. However, on time scales which are of the order of 
the energy current relaxation time some wave-like tran- 
sient behavior could be expected. However, a rigorous 
derivation of a macroscopic transport law from the micro- 
scopic equations of motion, especially for small systems 
in a low-dimensional environment, remains a formidable 
challenge. 

To conclude this Section we mention that the we limit 
ourselves here to the case of a single-well on-site poten- 
tial. The double-well case has been studied in Ref. 
A detailed analysis of second sound propagation in the 
three-dimensional lattice is given in Ref. llTI . 
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FIG. 1: Power spectrum of energy current for a chain of = 
4096 C = 1 and diflterent energies. Dashed lines are best fit 
with a Lorentzian line-shape a/(l -f r^oj^). 



III. LINEAR RESPONSE 



this means choosing a single-pole response of the form 



At a coarse-grained level, the most general linear con- 
stitutive relation between the energy current J(a;, t) and 
the local temperature gradient can be written as an ex- 
tension of the usual Fourier law (we limit ourselves for 
simplicity to a one-dimensional case). 



J{x,t) 
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(3) 

Substituting into the continuity equation, and assuming 
a decay of the kernel K at infinity, one obtains 



dT 
1)1 



^ = / dx' 



dt'K{x~x',t-t')^{x',t') (4) 

-or, OX'^ 



(up to the heat capacity). The generalized heat diffu- 
sion equation is formally solved by means of the Fourier- 
Laplace transform, with the convention 



yielding 



1 



dx / dt/(a;,i)e-*(«^-^*) , (5) 

oo Jo 



T{q,z) = 



To{q) 



—iz + q^K{q, z) 



(6) 



with To being the initial condition. The standard dif- 
fusive pole is recovered for a constant K(q, z) = D, i.e. 
for an instantaneous response, with D being the thermal 
diffusivity. The simplest improvement would be to in- 
clude a single relaxation time r, and no spatial memory 
effects, namely an exponentially decaying kernel in time 
of the Cattaneo-Vernotte type In the above notation 



Kiq,z) 



-iz + 1/t 



(7) 



The quantity v has the physical dimensions of a velocity, 
is defined as — 



(8) 



Substituting approximation ([7]) into Eq. ^ one can 
readily recognize that the resulting T{q, z) is nothing but 
the Laplace transform of the equation 



1 dT 
T~dt 



, d^T 
dx^ 



(9) 



This is known as the telegraph equation ^ . The equation 
contains an additional term with respect to the standard 
heat equation that affects the solution for times of order 
r. Indeed, the second-order derivative with respect to 
time leads to a finite velocity of perturbations. For in- 
stance, for a pulse initially localized at the origin, T = 
at distances > vt. For i 3> r, ordinary diffusive be- 
havior is recoveredi^. It should be also remembered, that 
Eq. ^ arises as a continuum limit of the persistent ran- 
dom walk (see Ref. [l^ and references within) . 

According to the fluctuation-dissipation theorem — 
the imaginary part of the response function K{q,Lo) is 
proportional to uj times the equilibrium power spectrum 
of the observable that couples to the external field. In 
the case of a temperature gradient, the observable to con- 
sider is the energy current, whose microscopic expression 
for model ([T]) reads — 
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FIG. 2; The relaxation time r as a function of energy for 
N = 2048, C = 1. The dashed line is a power-law fit. 



In order to test for the validity of approximation ([7]) 
we have first of all evaluated the power spectra of the 
q = component of the flux, j{q = 0,t) — X^nJ" from 
molecular dynamics. To this aim, we have performed 
microcanonical simulations by integrating Eqs. ([T]) (with 
periodic boundary conditions, Un = Uu+n) by means of a 
fourth-order symplectic algorithm ~. Initial conditions 
were chosen with the particles at equilibrium. Their ve- 
locities were drawn at random from a Gaussian distribu- 
tion and rescaled by suitable factors to assign the total 
energy per particle e to the prescribed value and to set 
the total initial momentum equal to zero. A suitable 
transient is elapsed before data acquisition for statisti- 
cal averaging. Conservation of energy and momentum 
was monitored during each run. The chosen time-steps 
(0.05-0.1) ensures energy conservation better then a few 
parts per million. The reliability of the spectra has been 
checked against different choices of the run duration and 
sampling times. 

Fig.[T]reports {\j{q — 0,w)p), where the average is per- 
formed over about 100 trajectories. The data show that, 
in agreement with Eq. ([7]), the spectra are fitted very well 
by a single Lorentzian over a wide range of energies and 
frequencies. The relaxation time decreases as a power of 
the energy density, r ^ with /3 « 1.4 (Fig. [2]). To 
compute the velocity v we measured independently the 
thermal conductivity k, by the standard Green-Kubo for- 
mula i.e. by integrating the flux autocorrelation in the 
time domain^. It turns out that the diffusivity D is ap- 
proximately equal to the thermal conductivity since the 
heat capacity is very close to one in our units (it varies 
only by a few percent in the considered energy range). 
As shown again in Fig. [2l k decreases as a power of the 
energy density. This is in agreement with previous work 
on related models — iS^. The corresponding exponent is 
very close to /3 meaning the that v as defined by Eq. ^ 



is roughly constant. 

In a second series of simulations, we computed the dy- 
namical structure factor, namely the square modulus of 
temporal Fourier transform of the energy density on the 
lattice 

en = Y + Y T 2"^ ~ ' ^ ' 



e{q,t) = ^ e„exp(-ign) , (12) 

n 

which is defined as 

Siq,u;) = {\eiq,Luf) . (13) 

The square brackets denote an average over a set of in- 
dependent molecular-dynamics runs. By virtue of the 
periodic boundaries, the allowed values of the wavenum- 
ber q are integer multiples of 2n/N . 

The data in Fig. [3] are representative of the numerical 
results. At low enough temperatures we see a peak at 
finite-frequency which suggests some kind of oscillating 
response i.e. the propagation of damped temperature- 
waves. Upon increasing e and/or decreasing q, the spec- 
tra display a central peak, akin to the one of an over- 
damped oscillator which signals the onset of diffusive be- 
havior. 

For not too large q the spectra are well fitted by a 
Lorentzian shape 



S{q,uj) 



So 



{u;^-Lui{q)r + {u;/riq)y 



(14) 



We found that the dependence of the parameters wo(<?) 
and T{q) from the wavenumber q is 



^0(17) 



r{q) 



1 



(15) 



where c and go are fitting parameters (see Fig.|4|). 

The observed form of the line-shape, Eq. (fT4|) . is fully 
consistent with what expected from Eq. However, 
the origin of the dependence of T{q) cannot be accounted 
for by the simplified kernel ([7]) and is presumably a sig- 
nature of spatial correlations. Nonetheless, we note that 
the value r obtained by extrapolating at g — > in the sec- 
ond of Eqs. is in very good agreement with the value 
measured from flux spectra (see again Fig ^ . Another 
important quantitative difference is in the characteristic 
velocities. In Fig. HJs we compare the velocity c from 
the fitting (|15p with v as given by definition ([8]). Since 
we found above that D and r are both proportional to 
(roughly) the same power-law , in Fig. |3] we report 
also the ratio of the corresponding proportionality con- 
stants. The measured values of v in Fig. ([4]) show that 
V < c in the considered range. 

We may surmise that the velocity c should be related to 
the group velocity of the harmonic waves. To account for 
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FIG. 3: Power spectrum of energy fluctuations for a 
N = 2048 C = 1 and (a) e = 0.1 (b) e = 0.5 and (c) 
Dashed lines are best fit with Ea. (|14|l . 
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FIG. 4: (a) Dependence of r(g) on the wavenumber for two 
different energy densities. Lines are best fit with the law pSI) . 
Triangles are the data for r reported Fig. [2] (b) Comparison 
of the relevant velocities emerging from linear response. The 
velocity c defined from Eq. (|14|) (dots) is very close to the 
maximum group velocity computed from the effective disper- 
sion, Eq. (llOp (solid line). The thermal velocity v is defined 
by Eq. ((Sjl (see text). Data are obtained for a chain of length 
iV = 2048 with C = 1. 



finite-temperature effects we have computed the renor- 
malized dispersion relation in the self-consistent phonon 
approximation — , 



where 



n'^iq) = w2(T)H-2C(l-cos(j) 



\T) = 1 



in') 



(16) 



(17) 



The second and fourth momenta of the displacements 
have been evaluated numerically. The data for c are very 
close to the maximum group velocity Cmax, namely the 
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largest value of <Kl{q)/dq as computed from formula ((T6)l . 
see Fig. g] 

We conclude that, although the telegraph equation Q 
accounts for the line-shape of the energy correlators, 
there are, at least in the considered times and length 
ranges, some quantitative deviations. The fact that v 
and c are different can be partly understood by noting 
that V is associated with diffusive processes and results 
from interaction of all possible lattice waves. The max- 
imum group velocity c will thus be an upper bound to 
V, but the two need not be equal^. It is interesting to 
mention that differences in the measured velocities were 
previously reported also in a model of hard-point parti- 
cles ^ . 



IV. HEAT PULSE PROPAGATION 

The linear response results reported so far suggest that 
the propagation of heat waves could be described at the 
macroscopic level, at least for small enough energies, 
through the telegraph equation ([9]). In order to test this 
conjecture, we performed nonequilibrium numerical sim- 
ulations where a heat pulse of temperature Tg is excited 
in a small region of width A7V within a chain otherwise 
at equilibrium at the bulk temperature < Tq and ob- 
served as it evolves. A typical pulse experiment proceeds 
as follows. First the chain is let evolve at a fixed en- 
ergy e — Th for a long enough equilibration time. Subse- 
quently, the portion [{N - AiV)/2, (TV + AN)/2] is put 
in contact with a Langevin thermostat at the tempera- 
ture To, while the rest of the chain is kept frozen in its 
equilibrium configuration. When the heated portion has 
reached thermal equilibrium with the thermostat, the lat- 
ter is removed and the equations of motion of the whole 
chain are integrated at constant energy. For the pulse 
propagation experiments we used the symplectic "Posi- 
tion Extended Forest-Ruth Like" (PEFRL) algorithm of 
Omelyan et al.-"* and a velocity Verlet algorithm for the 
microcanonical and Langevin integrations, respectively. 
Furthermore, since we used free-ends boundary condi- 
tions, wc took care to chose the integration time so as to 
avoid the propagating pulse to be reflected at the chain 
edges. In order to obtain highly accurate results, we av- 
eraged the time-dependent temperature profiles over a 
large ensemble of J\f independent realizations of the equi- 
librium configuration of both the bulk and the excited 
regions of the chain. All results reported in the following 
were obtained with Af = 2 x 10* and AN = 90, except 
where explicitly indicated otherwise. 

Typical spatio-temporal portraits of the temperature 
profiles describing the evolution of a heat pulse are re- 
ported in Fig. [5] for two different values of the bulk tem- 
perature. For small values of T},, the heat excitation gives 
rise to two symmetrical pulses, traveling in opposite di- 
rections at a constant speed. Such speed is extremely 
close to the maximum group velocity of the linear modes, 
Cmax = 0.618. As the bulk temperature is increased, the 



correlation time r is expected to decrease thus reducing 
the window of ballistic propagation. In fact, a diffusion- 
like evolution of the pulse is observed already at Tb = 0.04 
(Fig. [SId). However, it can be clearly appreciated that 
the bulk temperature is still not sufficiently large for the 
pulse to spread homogeneously at all heights. At a closer 
inspection, it is not difficult to realize that the the front 
baseline still clearly moves at a constant speed close to 

Cmax • 

Stated more precisely, our conjecture implies that, if 
the telegraph equation holds, then the motion of the heat 
pulse front should carry the information on the micro- 
scopic parameters that enter the continuum description, 
namely the speed v and the correlation time r (see again 
equation Let us imagine to cut the time-dependent 
temperature profiles that describe the chain dynamics af- 
ter excitation of the heat pulse at a given temperature T* , 
with Tfc < r* <IC Tq. Then, one may extract useful infor- 
mation by tracking the sites occupied by the pulse front 
at the height T* during its evolution. In other words, 
one may study the trajectories defined implicitly as 



T{x{t),t) = T* 



(18) 



where T{n,t) — {u^{t)) is the temperature field , the av- 
erage being computed over an ensemble of M independent 
trajectories. If T* is sufficiently close to the bulk tem- 
perature Th, the linear response results should hold, and 
the persistent random walk should be recovered, namely 



xit) 




(1 



-ft] 



(19) 



The analysis summarized in Fig. [6] proves the validity 
of our inference. Close to the background temperature, 
the propagation crosses-over from ballistic to diffusive 
on a time scale that is well predicted by the equilib- 
rium simulations (see Fig. [2]). For Tf, = 0.01, the equi- 
librium prediction would be a time scale of the order 
1/7 w 2 X 10*. Correspondingly, on our observation win- 
dow we observe purely ballistic propagation (upper right 
panel of Fig. [6]). Conversely, for Ti, = 0.3 the equilibrium 
relaxation time is about 150. Accordingly, we indeed ob- 
serve a crossover to diffusion within the time span of our 
simulations. From the fit we get 1/7= 127, in good 
agreement with the equilibrium prediction (lower right 
panel of Fig. [S]). At intermediate temperatures, an in- 
crease of the cut temperature height in the vicinity of 
Tfc causes the ballistic/diffusive cross-over to occur on 
shorter time scales, as expected since the relaxation time 
approaches the observation time. In fact, this can be seen 
as the very definition of the intermediate time scale. This 
is clearly illustrated in the middle right panel of Fig. [6] 
(Tfc — 0.1), where the transition becomes visible in the 
simulation time window upon raising the cut tempera- 
ture from T* = 0.104 to T* = 0.15. 

As a further check of consistency, we have examined 
the time-variations of the temperature field second mo- 
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FIG. 5: Spatio-temporal density maps of the normalized temperature profiles T{n,t) /To of two heat pulses propagating in 
chains at temperatures Tt = 0.001 (a) and Tt = 0.04 (b). Other parameters are: To = 20 x Tb, C = 1, N = 1500. 



ment /i2(i)j defined as 



(20) 



where z„ = w^/ ^ 



m^m and 



1. However, 

straight calculation of the second moment from its def- 
inition requires extremely accurate averages, due to the 
strong amplification of the fluctuations in the bulk re- 
gions away from (n) . For this reason, in order to monitor 
fi2 over time, two copies of the system were evolved start- 
ing from identical initial conditions but for the small re- 
gion around the center of the chain, where the heat pulse 
is generated. An accurate estimate of the pulse spread 
for finite values of Tf, can then be obtained from the sec- 
ond moment of the difference between the profiles of the 
two system clones. 

The result of such procedure is shown in Fig. [7] for two 
values of the bulk temperature in the case of a narrow 
excitation. As it shows, the heat pulse spreads according 
to the prescription of the telegraph equation, that is 



=M2(0) + — 
7 



7 



-7* 



(21) 



The time scales extracted from the fits are r — 917.5 
{Tb = 0.1) and T = 53.9 {Tb = 0.5) . These figures are 
in good agreement with the equilibrium relaxation times 
that we found by fitting the results of microcanonical 
simulations (data shown in Fig. [2]), t — 760.8 {Tb = 0.1) 
and T = 76.2 {Tb = 0.5). 

The best-fit values of the velocities v prove rather in- 
sensitive to the bulk temperature. We find v = 0.095 
{Tb = 0.1) and v = 0.1 {Tb = 0.5). Interestingly, these 
values are smaller than all estimates shown in Fig. 2] 
This is likely to reflect the temperature-dependence of 
the relaxation times at different heights within the pulse. 



For a finite temperature disturbance, the fronts do not 
spread at the same rate, causing the hotter portions to 
lag behind the advancing baseline. Overall, this should 
refiect in a lower value of the pulse velocity during its 
first ballistic stage. 



V. DISCUSSION 

In this report we have investigated the relaxation of 
temperature fluctuations in a discrete nonlinear chain. 
At low temperatures the relaxation time r of the energy 
current must be taken into account, leading to correc- 
tion to the standard diffusive behavior. Starting from 
numerical calculation of the response function and to the 
simplest level of approximation, we obtained the tele- 
graph equation, Eq. ([5]), and estimated the temperature 
dependence of its parameters. 

The comparison between the linear-response prediction 
and the noncquilibrium simulations reveals that the tele- 
graph equation provides a reliable macroscopic interpre- 
tative framework for quantifying the front propagation. 
In particular, we have shown that the second moment 
of the temperature field displays the ballistic/diffusive 
cross-over on a time scale in accordance with the correla- 
tion times of energy fluctuations extracted from equilib- 
rium simulations. The same is true for the front propaga- 
tion at temperatures close to the background. It might be 
surmised that this type of macroscopic description should 
apply to all models displaying normal energy transport, 
such as chains of coupled rotors or other one-dimensional 
models with pinning potentials^, or three-dimensional 
systems^^. 

By following the propagation at higher temperatures, 
the front appears to smear out in the course of time. It is 
likely that this effect could be captured by a kernel of the 
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FIG. 6: Analysis of heat pulse propagation. Spatio-temporal contour maps of the normalized temperature profiles T{n,t) /To 
(left panels) and time variation of the front baseline position at a given temperature height: T* = 0.0105 (b), T* = 0.104 and 
T* = 0.15 (d), T* = 0.35 (f). Parameters are. C = 1, U = 0.01 (a) and (b), Tt = 0.1 (c) and (d), Tt = 0.3 (e) and (f), 
To = 20 X Tb, C = 1, N = 1500. The straight lines in panels (b) and (c) are linear fits, while the other solid lines are fits with 
formula (IT9)) . 
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FIG. 7; Time evolution of the second moment of the temper- 
ature field for two values of the bulk temperature, Tt = 0.1 
(To = 0.5) and Tt = 0.5 (To = 1.0) (symbols) and two- 
parameter fits with expression (|21|l . The initial conditions 
/i2(0) have been fixed at the corresponding values extracted 
from the numerics. Other parameters are: C = 1, N = 5000, 
AiV = 10. 



Jeffreys type, as shown by Maruyama and co-workersi^. 
Another possible improvement would be to include spa- 
tial memory effects, by allowing for a space-dependence 
memory dependence in the kernel K, Eq.([3]). This is es- 



pecially important if one wishes to describe systems with 
anomalous transport properties^. Indeed, in this case 
heat propagation is quantitatively described by a Levy 
walk process^i^, which is precisely the generalization of 
the persistent random walk for of a memory decaying as 
a power law. As a consequence, the macroscopic equa- 
tion generalizing the telegraph equation should involve 
mixed spatio-temporal fractional derivatives Possi- 
ble nonlinear heat-wave propagation could also be taken 
into account by including higher-order powers of the lo- 
cal gradient. Further work along these lines is in progress. 

Finally, we remark that our results confirm that wave- 
like transport of energy is associated with the presence 
of optical branches in the linear spectrum. This could 
be of interest for heat transport in single-walled carbon 
nanotubes, where a significant contribution of optical 
phonons to wave-like conduction has been reportedi^. 
In this respect one may conjecture that our simplified 
model may provide an effective description of more com- 
plex quasi-one-dimensional nanostructures. 
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